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Abstract. We investigate the error implied by the use of the Zel’dovich approximation to set up the initial 
condit ions at a finite red shift Zi in numerical simulations. Using a steepest-descent method developed in a previous 
work (Valageas (2002a)) we derive the probability distribution V^Sr) of the density contrast in the quasi-linear 
regime. This also provides its dependence on the redshift Zi at which the simulation is started. Thus, we find that 
the discrepancy with the exact pdf (defined by the limit Zi —> oo) is negligible after the scale factor has grown by 
a factor a/ai ^ 5, for scales which were initially within the linear regime with Ci 0.1. This shows that the use 
of the Zel’dovich approximation to implement the initial conditions is sufficient for practical purposes since these 
are not very severe constraints. 
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1. Introduction 


In usual cosmological scenarios, large-scale structures 
in the universe form through the growth of small ini- 
tial density flu ctuations by gravitational instability (e.g., 
Peebles (1980) ). Besides, in most cases of cosmological 
interest the amplitude of these perturbations increases 
at small scales, as in the CDM model ( Peebles (1982)| ). 
This leads to a hierarchical scenario of structure forma¬ 
tion where smaller scales become non-linear first. They 
build small virialized objects which later become part of 
increasingly large structures. Thus, this halos give rise to 
galaxies or clusters of galaxies (depending on cooling pro¬ 
cesses). Unfortunately, this non-linear regime is very dif¬ 
ficult to handle analytically so that N-body simulations 
are a key tool to understand the formation of large-scale 
structures. They are even more important when one tries 
to follow the evolution of baryons which involves many 
processes (star formation, cooling, radiative ionization,..). 
Therefore, it is important to get a good estimate of the 
accuracy of such numerical simulations. 

In this respect, an obvious source of error is the gener¬ 
ation of initial conditions. Indeed, numerical simulations 
are initialized at a finite redshift Zi while they should be 
started at time tj = 0 when the relevant scales are exactly 
“linear”. In practice, one uses the Zel’dovich approxima¬ 
tion ( ZePdovich (1970)| ) to set up the initial conditions at 
Zj. This correctly matches the exact density and velocity 
fields at linear order but the higher-order terms are not 
exact (as compared with the fields defined by the same 
linear growing mode initialized at t = 0). This leads to 
a small error which can be made negligible by starting 


the simulations at a sufficiently large redshift Zj, when the 
scales of interest are far within the linear regime (e.g., 
puszkiewicz et al. (1995) , Baugh et al. (199^ ). 


These transients induced by the Zel’dovich approxi¬ 
mation were investigated in ^coccimarro (199^ for the 
first few order moments {q < 8) of the density and ve¬ 
locity helds using a perturbative approach. However, this 
approach involves lengthy calculations which worsen at 
higher orders and it is not obvious how to get an esti¬ 
mate of the error for the probability distribution func¬ 
tion (pdf) V{5r) of the density contrast from these mo¬ 
ments. Therefore, in this article we show how one can 
apply to this problem a steepest-descent method devel¬ 
oped in a previous work (Valageas (2002a)) which pro¬ 
vides a rigorous derivation of V^Sr) in the quasi-linear 
limit. This approach is quit e general since it applies to 
Gaussian initial conditions ( Valageas (2002a)) as well as 
to non-Gaussian initial conditions (e.g., Valageas (2002b) ) 
or to th e tails of P(Sr) in both linear and non-linear 
regimes ( Valageas (2002c) ). Also, it is actually quite in¬ 
tuitive. Here we describe how it can be used to derive the 
dependence on the initial redshift Zi of the error induced 
on the pdf 'P{6r). It also yields the moments of the density 
field at any order. 


This article is organized as follows. First, in Sect.||we 
recall the path-integral formulation which allows us to de¬ 
rive the pdf V{5r) in terms of initial conditions. Then, 
in Sect.^l we apply this approach to the gravitational dy¬ 
namics defined by Zel’dovich initial conditions at a finite 
redshift Zj. Finally, in Sect.^we present our numerical re¬ 
sults. 
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2. Steepest-descent method 

We first recall in this section the basis of the steepest- 
descent method built in Valageas (2002a)| . This also allows 
us to introduce our notation. We are interested in the 
statistics of the density contrast Sr within a spherical cell 
of physical radius R, volume V: 




( 1 ) 


AL(ri,r2) = (5L(ri)5L(r2)). 


( 2 ) 


= J 


dk PL{k). 


(3) 


dri dr2 


(4) 


Sion: 


V’(y) = 

9=0 


tvTiM 

q\ cr 29 ■ 


( 6 ) 


= with = 

9=2 


since we have (5 _r)c = 0. 

Then, the Gaussian average which appears in eq.(|^ 
can be expressed as the path-integral: 


i,{y) = (DetA^i) ^ j [d5L(r)] 
where we introduced the action S'[5 l]: 


= y Sr[5l] 


\R) 


Sl.AT^.Sl 


(9) 


( 10 ) 


Note that in this article all length-scales will be physical 
lengths (and not comoving coordinates). Here ^(r) is the 
non-linear density contrast at the physical coordinate r, 
at the time of interest. It is fully determined by the linear 
growing mode 5L(r) since we assume that the decaying 
mode has had time to vanish. Moreover, for Gaussian ini¬ 
tial conditions the statistics of the linear density field ^L(r) 
are fully defined by the kernel Ai(ri,r 2 ): 


Note that AL(ri,r 2 ) = Ci(|i "2 — ri|) is simply the two- 
point correlation function of the linear density field, which 
is related to the linear power-spectrum PL(k) by: 


Here is the non-linear density contrast 5r within 

the spherical cell V obtained for a linear density field 
5l{^)- The action S[6l] is independent of the normaliza¬ 
tion of the power-spectrum PL^k) since A^ oc cr^. Then, 
it is clear that the path-integral in eq. (|) is dominated 
by the minimum of the action S in the limit a 0 for 
a fixed y. Indeed, the contributions from other states Sr 
are exponentially damped relative to this point. Moreover, 
the steepest-descent approximation becomes exact in this 
limit. This allows us to compute the generating functions 
ipiy) and 'tjj{y) in this limit. Finally, we obtain the pdf 
V{dR) from eq. through the inverse Laplace transform: 


/ +ioo 

-ioo 


dy 


2'!ria‘^{R) 


,lySR-ip(v)]/P(R) 


( 11 ) 


Besides, the rms linear density fluctuation a{R) in a spher¬ 
ical cell of radius R, volume V, is: 


a\R)^{6lR)= 

where Sr^r is the linear density contrast within the spher¬ 
ical cell of radius R. 

Next, in order to derive the moments and the pdf 
of the density contrast it is convenient to introduce 
the generating function tpiv) defined by the relation: 

/ OO 

ddfl V{Sr). (5) 

The average (..) in eq.(||) denotes the mean over the 
Gaussian initial conditions, which are defined by the lin¬ 
ear growing mode Sr^t). The last equality defines the 
pdf V{Sr) of the density contrast. Thus, tp{y) is also the 
“rescaled” Laplace transform of P{SR).lt also yields the 
moments of the density contrast Sr through the expan- 


Next, in order to apply the steepest-descent method 
we need to derive the saddle-point of the action As 

shown in Valageas (2002a) a great simplification occurs 
thanks to the spherical symmetry of the problem. This 
symmetry arises because the initial conditions (i.e. the 
linear growing mode Sr) are statistically homogeneous and 
isotropic and we consider the density contrast Sr within 
spherical cells. Then, one obtains a spherically symmetric 
saddle-point. Besides, for spherically symmetric states we 
can write: 


Sr = ^[Sl,Rl] 
i?i = (1 + Sr)R^ 


( 12 ) 


where the function P[Sr^Rj^] is given by the usual spher- 
ical collapse solution of the equations of motion, see 


Peebles (1980)| and [Valageas (2002a)l The second equa¬ 


tion in w merely expresses the co nservation of mass. 
Then, as shown in Valageas (2002a) from eq.(|l^, for a 
power-law linear power-spectrum: 


PR{k) oc fc" with — 3 < n < 1, 


(13) 


the generating function (p{y) given by the steepest-descent 
method is defined by the implicit system: 


It is convenient to introduce also the generating function 
of the cumulants defined by: 

i>(y) = e-^d!)lP ^ 

which exhibits the expansion: 


T = -y G'{t) 


(14) 


(7) 

[ v(?/) = y Q{t) + — 



Here we introduced the function Q{t) which is 
the function J-[Sr^r^] defined in eq.(|T^) by: 

related to 

(8) 

e;(r) =^[-r (l + g[r])-("+3)/6 . 

(15) 
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This yields the pdf V{5r) through eq. ( 0 ) and eq. ( 0 )- 
We refer the reader to Valageas (2002a)| for a detailed dis¬ 
cussion of the steepest-descent method within this frame¬ 
work. 


3. Zel’dovich initial conditions 

The results described in the previous section were de¬ 
rived for Gaussian initial conditions at the time U = 0. 
Indeed, eq. ( 0 ) implies that Sr/6l,Rl 1 for t ^ 0. 
However, in numerical simulations the initial conditions 
are set up at a non-zero time ti > 0. This leads to a 
small error for the properties of the density held. This 
discrepancy vanishes in the limit ti —> 0. Here we are in¬ 
terested in the quantitative estimate of this error for the 
moments and the pdf V^Sr) of the density contrast. This 
can easily be done through the steepest-descent method 
recalled in SectJ^Indeed, all the steps involved in the 
derivation of eq.(p^ and eq.(|l^ remain valid because the 
initial conditions are again statistically homogeneous and 
isotropic and they are still dehned by the Gaussian lin¬ 
ear growing mode Sl{y). However, the function J^[5l,Rl] 
which describes the spherical dynamics is no longer given 
by the usual spherical collapse equations (i.e. eq.(39) and 
eq.(40) in Valageas (2002a)). Therefore, in order to derive 
the pdf 'P{5r) we simply need to obtain the new function 
^[5l,Rl\- 


3.1. Characteristic function iiF[6L,Ri^] 


In the following we consider a critical-density universe 
where the analytic expression of the spherical dynamics 
is quite simple. However, our method obviously applies 
to any cosmological parameters. One simply needs to use 
the relevant expressions for the (time-dependent) func¬ 
tion iF[SL,RL]- Let us first consider the case of a spherical 
overdensity Sr^r,^ > 0. Thus, we study the evolution of a 
spherical region of constant mass M and physical linear 
radius RL{t) oc a(t), where we note a{t) the scale-factor. 
We can take: 

RL{t)=a{t)=t^^^. (16) 


The use of the Zel’dovich approximation to set up the 
initial conditions at ti > 0 does not break the spherical 
symmetry so that the equation of motion is given by the 
spherical dynamics: 


_ 


(17) 


Here we introduced Newton’s constant t/N and the actual 
physical radius R{t) of the spherical overdensity. For a 
positive overdensity the general solution of eq. 0) reads: 

{ R = ^(1 ~ cos 9) 

< with = O^MB^. (18) 

t = C + B(9 — sin0) 

Note that the offset C in eq.(p^ vanishes in the limit 
ti —> 0 where we must recover the usual spherical dynam¬ 
ics. However, it must be taken into account in order to 


write the general solution to eq. ( p^ . In particular, the er¬ 
ror introduced by the Zel’dovich approximation at ti > 0 
yields a small non-zero value for this constant of integra¬ 
tion C. 

Next, in order to obtain the function TISr^Rj^] we need 
to derive the parameters H, B and C. First, we note that 
for a critical-density universe we have: 

^ = irk?' ' 9 I 7 

where we used eq. Thus, we can express the param¬ 
eter A in eq. in terms of B. Next, we must specify 
the initial conditions at time f. We consider an overden¬ 
sity with a linear density contrast at the time t which 
we are interested in. At time ti the linear density con¬ 
trast was dR^i = {ai/a)SR (where we note Oi = a{ti) and 
a = a{t) the scale-factors at times ti and t). The Zel’dovich 
approximation means that we apply to the particles the 
displacement given by linear theory. The exact radius R 
obeys R = Rr{1 + Sr)~^^^ therefore the Zel’dovich ap¬ 
proximation yields: 

i?i = Oi ^1 — 

where we used the normalization 0)- From eq.(p0|) we 
see that the physical velocity v{t) = dR/dt at time ti is: 



Vi = Hitti 



( 21 ) 


within the Zel’dovich approximation, where we used 
SR(t) oc a(t) and we introduced the Hubble constant 
Hi = H{ti) at time ti. Since a{f) = we also have: 
H{f) = aja = 2l{3t). The initial conditions ( |^ ) and (21) 
at time f define the constants of integration B and C in 
eq.(p^) through: 

/2 

Ri=l-B'^j (l-cos6»i) 


< ti = C + B{9i-sin9i) (22) 

/ 2 sin0i 

\ 9^ / 1 — cos 9i 


This system also determines the parameter 9i associated 
with the time ti. In particular, from eq.(^) we see that 
the combination Rivf only depends on the angle 9i. Using 
eq.(0) and eq.(0) we obtain: 


Rivf = -(l-|-cos0i) 




(23) 


where we used SR^i = {ai/a)SR. This yields 0i as a func¬ 
tion oi Sr. Substituting into the first line of eq.(0), using 
eq-d^), we get: 


B = 





1 Oi 

3 a 


Sr 


3/2 


X 


1 ^ -L ^ 

1 “ 77 — ol + 77 
15 a 45 



(24) 
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Next, the angle 0 which corresponds to the time t of in- the Zel’dovich approximation. It is directly obtained from 

eq.(|^) and the first equality in eq.(| 2 ^) which gives (for 
any positive or negative Sl)- 


terest is given by: 

t — t[ = B[{9 — sin0) — {9i — sin0i)] 
while the density contrast Sn obeys: 


(25) 


1 + iJr = 


R 


= a 


3 / 2^2 

9 


(1 — cos 9) 


-3 


Putting together eq.(|23D, eq.(|25|) and eq.(| 
eq.(24), we finally obtain, for > 0: 


cos0i = l-^-<5L+ 

6 a \6 a ) \6 a 


sin 0 = 0 ; — sin 0 ; -|- 
- 10^/2 




1 - 
3/2 


© 


3/2' 


d a 


1/2 


-3/2 


( 1 ^ “i;: , 4 /a; y 

15 a 45 (a 


3/2 


1 + = 


3 (1-(1-COS0) 


cosh 77 ; = 1 - ^ — Sl 
c> a 




sinh rj — rj = sinh 77 ; — 77 ; -|- 

(t) 


1 - 




(7) 


3/2' 


0 


1/2 


1 - -^—Sl 

c> a 


-3/2 




-3 


- 1 


(29) 


(26) where the superscript “Z” refers to the Zel’dovich approx¬ 
imation. Next, we can obtain the coefficients defined 
using in eq.(||) by expanding the functions R^{Sl), (t) and 

ip^{y) around zero. This yields: 


(27) 


This system of equations fully defines the function 
T[5l,Rl\ ctf time t for positive overdensities since we have 
6 r = R{ 6 l)- Thus, the first line gives 0;, the second one 
yields 0 and the last one provides 6 ^. Note that for i5l ^0 
we have 0; —> 0 and 0 —> 0. Moreover, since we con¬ 
sider the quasi-linear regime we can restrict ourselves to 
6l < 3/2 and the angle 0; obeys 0 < 0; < tt while we have 
0 < 0 < 27r. Of course, we can proceed in the same fashion 
for underdensities which yields, for 5l < 0 : 


(28) 


1 ^ “‘A 

1 - TT —© 

15 a 


4 

45 


or 


3/2 


1 + = 


10 ' 






-I 3 


3 {^-^^^L)icoshr]-l) 


Here we also have 77 ; > 0 and 77 > 0. Finally, as explained 
above, the determination of the function J-[ 6 l,Rl\ also 
gives the pdf V{6r) through eq. (0), eq.dil) and eq. ©• 

3.2. Zel'dovich pdfV^{5R) 

At this point, it is interesting to consider the pdf V^{Sr) 
at the initial time t;. This is in fact the pdf defined by 


G^{r) = ip^iy) = -y-f 

which gives for the skewness: 

= 4 — (n -I- 3). 


"'■‘-^+..(30) 


(31) 


Of course, we recover the results obtained in 
[Bernardeau fc Kofman (1995) who derived the pa¬ 
rameters Sq in the quasi-linear limit defined by the 
Zel’dovich approximation using perturbative 
Indeed, as shown in Valageas (2002a) the 


defined by 

means, 
steepest- 


in Valageas (2002a) 
descent method recalled in Sect J| gives the same results 
for these coefficients Sq as the perturbative technique 
pr esented in Bernard eau (1992,1994). We refer the reader 
to Valageas (2002a) (Sect.4.1) for a detailed comparison 
of the steepest-descent method with the perturbative 
approach. 

Obviously, since the function R^{Sl) defined by the 
Zel’dovich approximation is different from the exact spher¬ 
ical function R{6 l), the pdf V{Sr) and the coefficients Sq 
obtained at time <; differ from the value defined by the 
exact gravitational dynamics (where the initial conditions 
would be set at t =^- For instance, the skewness factor 
S^ obtained in eq.(^) is different from the e xact value 
S 3 = 34/7 — (ti -I- 3) (e.g., Bernardeau (1994) ). However, 
as the initial conditions are set up at earlier times t; < 1 / f 
the error induced by the non-zero value of t; decreases 
and in the limit t[ 0 the skewness factor measured 
on linear scales at time t tends towards the exact value 
S 3 = 34/7 — (77-1-3). The goal of the present article is 
precisely to study how fast this relaxation proceeds. In 
other words, we wish to see at which redshift 2 ; we must 
set up the initial conditions so that the error due to the 
use of the Zel’dovich approximation at t\ can be safely ne¬ 
glected. This can be obtained fro m th e time-evolution of 
the function J-{Sl) derived in Sect.3.1. 


3.3. Lowest order moments 

The eq.(^) and eq.(^) provide a parametric definition of 
the function tF{6]3) at any time t > t;. Besides, in the limit 
^ 0 we have 0 ; ^ 0 , 0 —> 0 (or 77 ; —> 0 , 77 —> 0 ) and 
5r 0. Therefore, by expanding the relations in eq.(p^ 
and eq.(p^ around zero we obtain the Taylor expansion of 
p{y) around 1 / = 0. This yields the coefficients Sq through 
eq.(^). These expansions are straightforward with a sym¬ 
bolic mathematical package (e.g., Mathematica), up to the 
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required order. For instance, we give below the expansions 
up to the second-order which is needed for the skewness. 
We get: 


^{^l) = 5l + 

and: 

Q{t) = -t + 
which gives: 

‘fiy) = -y + 

and: 


17 1 Oi 2 ^ai\ 

^ ~ S'a ^ WJ 


7/2' 


13 

42 


n 

6 


1 Oi 

5 a 


2 


0 ) 


Si 


7/2' 


13 

42 


n 

6 


1 Oi 2 

57 ^ 


e) 


7/2' 


(32) 


(33) 


.. (34) 


C 34 . doi 

53 = y - (n -k 3) - 


(?) 


7/2 


(35) 
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We can chec k that our result d^ ) agrees with the value 
obtained by Scoccimarro (1998)| from a perturbative ap¬ 
proach. The advantage of our formulation is that it pro¬ 
vides at once the coefficients Sq up to any order. Besides, 
it also gives the pdf V{5r) through eq. ( 0 )- Of course, we 
can see that for Uj = a we recover the predictions of the 
Zel’dovich approximation discussed in Sect. |3.2| while for 
Oi —> 0 we recover the exact values of the gravitational 
dynamics derived in Valageas (2002a)[ 


Fig. 1. The dependence on Z{ of the ratio 53 (^ 1 ) 753 ( 00 ) 
of the skewness factor measured at z = 0 with respect 
to the exact value (obtained when Zi = 00 ). We consider 
three linear power-spectra (i.e. n = - 2,-1 and n = 0 , 
from top to bottom) in a critical-density universe. The 
coefficient is obtained from eq.(^^. The data point 
shows the measure of the skewness in a numerical simula¬ 


tion (Bernardeau (1994)) with zi = 8.7 for n = —1. 


results only hold on large scales with a ^ 0 . 8 , which be¬ 
long to the q uasi-linear regime. T hese results were already 
obtained in Scoccimarro (1998) from a perturbative ap¬ 
proach hence we shall not comment them further here. 


4. Numerical results 

Finally, we present in this section the numerical results 
obtained from the analysis described in SectJ|. We com¬ 
pare our results with the numerical simulations analysed 
in Bernardeau (1994)| . They use an adaptative P^M code 
with 2.1 X 10® particles. The error bars are estimated by 
dividing the simulation box in eight equal subsamples and 
by making eight different measurements. 


4.1. Skewness 

First, for completeness we display in Fig|| the dependence 
of the skewness factor S 3 measured today (i.e. at z = 0 ) 
with the redshift Zj at which the initial conditions are set 
up. We consider three linear power-spectra: n = -2,-1 
and n = 0. The curves show the ratio 53 (zi)/ 53 (c»), where 
53 (c») is also the exact value of the skewness factor on 
quasi-linear scales (i.e. when the initial conditions are set 
up at t = 0). Of course, we can check that the error over 
S 3 decreases as z, —> 00 . Besides, we can note that the 
discrepancy with the exact value grows for larger n. The 
data point shows the measure of S 3 in the numerical sim¬ 
ulation for n = — 1 after the scale-factor has grown by 
a factor a/oj = 9.7 (i.e. Zi = 8.7). We can see that the 
error due to the Zel’dovich approximation at the initial 
time is already negligible, as compared with the accuracy 
limitations due to finite size effects. 

As seen in pcoccimarro (199^ the error also increases 
for larger order q of the parameters Sq. Note that these 


4.2. Probability distribution V{5r) 
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Fig. 2. The pdf V{6r) for n = —1, = 1 and cr = 0.74, 

at redshift z = 0. The solid line shows the “exact” pdf 
(i.e. when the initial conditions are set up at Zi = 00 ). 
The dashed line corresponds to Zi = 5 and the dotted line 
to Zi = 0 (i.e. to the Zel’dovich approximation). The data 
points (obtained from numerical simulations) are taken 
from Bernardeau (1994) with Zi = 15. 


We can note that the lowest order moments do not 
fully describe the pdf V{5r) so that it is not easy to 
get an accurate estimate of the deviation of V{Sr) from 
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its asymptotic value from the first few coefficients Sq. 
Fortunately, as detailed in SectJ| we also derived the pdf 
T’{5ji) itself. Thus, we present in Fig.^the dependence on 
Zi of the pdf V{5}i) measured today at z = 0, in a critical- 
density universe with n = — 1. The case z, = 0 actually 
gives the reconstructed pdf implied by the Zel’dovich ap¬ 
proximation. We can clearly see the difference with the ac¬ 
tual pdf which comes from the exact gravitational dynam¬ 
ics with initial conditions set up at ti = 0 (obtained within 
the quasi-linear approximation). However, we note that at 
Zi = 5 the pdf is already very close to this asymptotic value 
and at Zi = 10 it would be nearly indistinguishable in the 
figure. Therefore, we can conclude that for practical pur¬ 
poses the scales which obey Ci ^ 0.1 at the time ti when 
the initial conditions are set in are well described by the 
numerical simulation after the scale-factor has grown by a 
factor a/oi ^ 5. In particular, it is clear from Fig.^ that, 
provided these constraints are satisfied, the error bars of 
the numerical simulations are largely dominated by other 
problems (e.g., finite size effects, finite resolution). Note 
that the measures obtained from the simulations corre¬ 
spond to Zi = 15 where the error due to the set-up of the 
initial conditions is indeed negligible. 


For simplicity we have only considered the case of a 
critical-density universe so far. Nevertheless, it is clear 
that our approach can be extended to any cosmology. One 
simply needs to use the relevant equations of motion which 
describe the spherical dynamics of the saddle-point (e.g., 
eq.(@). This merely changes the characteristic function 
However, this is not necessary for practical pur¬ 
poses. Indeed, since observations show that the present 
universe is close to the critical density (Hm ^ 0.2) at 
high redshift we have Hm(z) ~ 1 to a very good accu¬ 
racy. Moreover, in order to study large-scale structures in 
the present universe numerical simulations must start suf¬ 
ficiently early at z > 20, especially if one is interested in 
the behaviour of baryons which involves early processes 
like reionization at z ^ 6. Therefore, during the initial 
stages of realistic simulations we have indeed Hi„(z) ci 1. 


Nevertheless, in case one is only interested in the 
present universe at z ~ 0 it is important to take 
into account the effect of the cosmological parame¬ 
ters on the time-evolution of the transients due to the 
Zel’dovich approximation. To do so, one must use the rel¬ 
evant characteristic function T[5l,Ri,] which is no longer 
given by eq.(p7|) and eq.(|^). For instance, if Ha > 0 
an additional term appears in the equation of motion 
fi. which we must integrate. However, as noticed in 
Nusser & Colberg (1998) the evolution of density pertur¬ 
bations is almost independent of the cosmological param¬ 
eters if we use the linear growth factor D+{t) as a time 
variable. This means that the results of Figj^ remain valid 
except that the redshift Zi = 5 shown by the dotted line 
actually corresponds to the redshift z[ given by: 


^+( 0 ) 


1 + Z[, 


where we used the fact that for the critical-density uni¬ 
verse we have D+ oc a oc (1 -|- z)~^. For instance, for a 
low-density flat universe with Hm = 0.3 and Ha = 0.7 we 
get: 

= 5 ^ z( = 6.7 for Hm = 0.3, Ha = 0.7 (37) 

Thus, as expected we find that a larger number of expan¬ 
sion factors are required in the case of low-density uni¬ 
verses as structure formation is frozen at late time by the 
more rapid expansion of the universe. 

5. Conclusion 

Thus, in this article we have shown how to derive the pdf 
V{5r) of the density contrast implied by the use of the 
Zel’dovich approximation to set up the initial conditions 
at a finite redshift Zi. This allows us to obtain a rigorous 
quantitative estimate of the error due to this effect for 
the pdf V{5r) and the moments of the density contrast. 
Our results apply to large scales which are still in the 
quasi-linear regime. This is quite sufficient for practical 
purposes since it is clear that it makes no sense to use the 
Zel’dovich approximation at time ti for scales which are 
already non-linear. 

Then, we found that the error for the pdf V{5b) associ¬ 
ated with the Zel’dovich approximation at z, is negligible 
after the scale factor has grown by a factor a/oi ^ 5, 
for scales which were initially within the linear regime 
with (Tj ^0.1. In low-density universes one needs a larger 
expansion factor ratio a/aj. In fact, to a good approxi¬ 
mation the relevant quantity is the growth factor ratio 
Zl+(t)/Zl+(ti), which should be larger than 5. Thus, in or¬ 
der to use numerical simulations it is important to check 
that these constraints are indeed satisfied over the scales 
one is interested in. However, these conditions are not 
very severe, hence it appears that this effect is not the 
main source of error. Therefore, it is sufficient to use the 
Zel’dovich approximation to set up initial conditions and 
one does not need to use higher-order approximations. 
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